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Abstract 

We present lattice results for spin-1/2 fermions at unitarity, where the effective range of the 
interaction is zero and the scattering length is infinite. We measure the spatial coherence of 
difermion pairs for a system of 6, 10, 14, 18, 22, 26 particles with equal numbers of up and down 
spins in a periodic cube. Using Euclidean time projection, we analyze ground state properties and 
transient behavior due to low-energy excitations. At asymptotically large values of t we see long- 
range order consistent with spontaneously broken U{1) fermion-number symmetry and a superfluid 
ground state. At intermediate times we see exponential decay in the t-dependent signal due to 
an unknown low-energy excitation. We probe this low-energy excitation further by calculating 
two-particle correlation functions. We find that the excitation has the properties of a chain of 
particles extending across the periodic lattice. 
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I. INTRODUCTION 

The unitary limit describes a many-body system of nonrelativistic spin-1/2 fermions 

with zero-range attraction and infinite scattering length. While the unitary limit has a 
well-defined continuum limit and strong interactions, at zero temperature it has no intrinsic 
physical scale other than the interparticle spacing. This implies for example at zero tem- 
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perature the energy per particle and pairing gap are both given by the Fermi energy times 
a dimensionless constant. 

The universal nature of the unitary limit endows it relevance to several areas of physics 
and the subject has received much recent interest. The ground state is believed to be 
superfluid with properties somewhere between a Bardeen-Cooper-Schrieffer (BCS) fermionic 
superfluid at weak coupling and a Bose-Einstein condensate (BEC) of bound dimers at strong 

n n n 

coupling pj, 121, |3|. In solid state physics it has been suggested that the crossover from 
BCS fermionic superfluid to BEC bosonic superfluid describes the pseudogap phase in high- 
temperature cuprate superconductors |4|. In atomic physics BCS-BEC crossover has been 
studied extensively with trapped ultracold ^Li and atoms. The atoms are sufficiently 
far apart that the effective range of the interaction is negligible while the scattering length 
can be adjusted using a magnetic- field Feshbach resonance ABB , 0] . In nuclear physics 
the unitary limit is relevant to the properties of cold dilute neutron matter. The neutron 
scattering length is about —18 fm while the effective range is 2.8 fm. Therefore the unitary 
limit is approximately realized when the interparticle spacing is about 10 fm, roughly 0.5% 
of normal nuclear matter density. Superfiuid neutrons at around this density may be present 
in the inner crust of neutron stars joifiol. 

In this study we measure the low-energy states of unpolarized spin- 1/2 fermions in the 



ll\ to measure the 



unitary limit. We use the same lattice projection technique used in 
ground state energy. We start with a quantum state with the desired quantum numbers 
and use the Euclidean time projection operator e"^* to filter out high-energy excitations. 
We work with finite systems where the energy spectrum is discrete. Therefore we don't 
have gapless modes which appear only in the thermodynamic limit. After this filtering we 
measure spatial correlations of the superfiuid order parameter as a function of the projec- 
tion time t. At asymptotically large values of t we see long-range order consistent with 
spontaneously broken U{1) fermion-number symmetry and a superfiuid ground state. At 
intermediate times we see exponential decay in the t-dependent signal due to an unknown 
low-energy excitation. We probe this low-energy excitation further by calculating two- 
particle correlation functions. We find that the excitation is consistent with a quasi- ID 
subsystem of particles extending across the periodic lattice. 
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II. OFF-DIAGONAL LONG-RANGE ORDER 



Written in continuum notation the Hamiltonian for the unitary hmit is 

H = ^ J £ral{r)V^a^{f) + C J £r a\{r)a\ir)a^ir)ai{r). (1) 

tti and a J are annihilation and creation operators for fermions with spin i. The mass of 
the fermion is m, and the coefficient C is cutoff dependent. We discuss later how in lattice 
regularization C is tuned to make the s-wave scattering length infinite. 

The unitary limit Hamiltonian has a global U{1) fermion-number symmetry 

(2) 



'at(r)" 













where is any real constant. It also has a global SU{2) spin symmetry 



"at(f)" 













(3) 



where a denotes the 2x2 Pauli spin matrices and is any real constant three-component 
vector. Since there is no coupling between intrinsic spin and orbital angular momentum, 
this SU{2) symmetry should be regarded as an internal symmetry decoupled from spatial 
rotations. 

The lowest-dimensional local bosonic operator that can be constructed from the annihi- 
lation field operators is 

^l)'^{r) = a^{r)ai{r). (4) 

We note that ip'^ is invariant under the SU{2) spin symmetry but phase rotates under the 
f/(l) fermion-number symmetry, 



(5) 



Therefore if there is some critical temperature below which ip"^ has long-range spatial corre- 
lations, 

hm U'\f)^\6)) ^ 0, (6) 

then the f/(l) fermion-number symmetry is spontaneously broken. While there are many 
different ways to characterize superfluid behavior, this condition of off-diagonal long-range 



order 



12 



13| is usually regarded as the standard definition for superfiuidity. 
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III. MEASURED OBSERVABLES IN CONTINUUM NOTATION 



In order to highlight the physics content of the lattice calculation we summarize the 
measured observables in continuum notation. We refer to a state with N-^ up-spin fermions 
and down-spin fermions as an N^, state. We also specify the total momentum P and 
total spin S of the SU{2) spin representation. Let j^E'o'"^) be the free Fermi ground state 
for the A^, system with P = and S = 0. We filter out high-energy states using the 
Euclidean-time projection operator e~^*, 

= e"^* K'^") . (7) 
With the projected states we measure the correlation function G^2{f,ti,t2), 

where F is a renormalization coefficient set by the condition 

lim G^2(0,ti,t2) = 1. (9) 

il ,i2— ►oo 

We consider only finite systems where the energy spectrum is discrete. For large t we 
have the asymptotic form 

\^>{t)) = coe-^«* l^o) + cie-^^* + ■ • • . (10) 

|\l/o) is the normalized ground state with energy Eq, and |\E'i) is the normalized first excited 
state with energy Ei. We anticipate the possibility of degeneracy in the low-energy spectrum 
corresponding with two-particle excitations with momenta k and —k. Since we use a periodic 
cube with lattice regularization, quantum states in the same irreducible representation of the 
cubic lattice symmetry group 5*0(3, Z) are exactly degenerate in energy. But there might 
also be approximately degenerate excited states whose energy differences are not adequately 
resolved by the finite values of t measured. We use the notation for the complete set 

of normalized degenerate excited states with energy Ei. Using the asymptotic form ffTOj) 
we get 

G^. (f, ti, h) = Aooif) + Aoi(f)e-(^-^»)*^ + Aio(f)e-(^^-^°)*^ 

+ An(r)e-(^i-^°)(*^+*^) + --- (11) 
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where 

Aoo(r1 = i(*o|V''^(r>'(0)|*o), (12) 

r V l^ol 

Ao(r) = ^ E V''^(r1V'(0) |*o) , (14) 

1 ^ \co\ 1 V |co| 

From the definition of F it follows that Aoo(O) = 1. In order to determine whether 
the strongest transient signal comes from Aoi{f) and Aio{f) or Aii{f) we will calculate 
G^2{f, ti, ^2) for two different large time combinations. In one case we let ti — t2 and send 
both to infinity. In the other case we let t2 be large but fixed and take the limit as ti goes 
to infinity. By analyzing the time dependence we extract the energy difference Ei — Eq, 
and from the spatial dependence we measure the momentum of the particles in which 
couple to '0^^(r)V'^(O). 

IV. LATTICE FORMALISM 

Throughout our discussion of the lattice calculation we use dimensionless parameters and 
operators which correspond with physical values multiplied by the appropriate power of the 
spatial lattice spacing a. When we need to specify quantities in physical units wc write the 
subscript 'phys'. In our notation the four-component integer vector n labels the lattice sites 
of a 3 + 1 dimensional lattice with dimensions x Lt. fis gives the spatial part of n, while 
rit is the time component. We write = 1, 2, 3 for the spatial lattice unit vectors and for 
the temporal lattice unit vector. The temporal lattice spacing is given by a^, and at — at/ a 
is the ratio of the temporal to spatial lattice spacing. We also define h — at /{2m), where 
m is the fcrmion mass in lattice units. 

We briefly discuss four different lattice formulations: the transfer matrix formalism with 
and without auxiliary fields, and the path integral formalism with and without auxiliary 
fields. While in the main calculation we use only the transfer matrix formalism with auxiliary 
fields, the other formulations are useful to provide numerical checks of the simulation data. 
The four formulations agree exactly even for nonzero spatial and temporal lattice spacings. 
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The first formulation we consider is tlie patli integral formalism with auxiliary fields. This 
has been used in varying forms in several grand canonical mean field calculations and lattice 
simulations at nonzero temperature [li 
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mm. We let cM) and c* n) be 



anticommuting Grassmann fields for spin i, and s{n) be an auxiliary Hubbard-Stratonovich 
field. Let Z be the partition function 



Z (X J DsDcDc* exp [—S{s,c, c*)] 



where 



Ds = Y\ ds{n), 

n 

DcDc* = Y[dci{n)dc*{n). 

n,i 

We take as our standard lattice path integral action 

S{s,c,c*) =J2 k(n)Q(n + 6) -e^^^^(")+^(l -6/i) 



(16) 

(17) 
(18) 



cAn cAn 



h [c* (n)Q(n + /,) + c*{n)ci{n - /,) + ^ t^^^^]^ • 



(19) 



nis.i 



If we include a chemical potential then the action becomes 



S{s, c, c*) = \c*{n)c,{n + 6) - eV^^«(")+^e^"*(l - Qh)c*{n)ci{n) 

n,i 

- /ie^"' Y k(r^)Q(n + I) + c*(n)Q(n - q] + ^ [sin)] 



(20) 



ft.ls 



From this point on it is easy to follow how the chemical potential enters into the other 
formulations. Therefore we simplify the discussion by setting the chemical potential to 
zero. 

In order to connect the path integral with the transfer matrix formalism, we use the 



correspondence |21 
Tr 



4(nJ,ai(ns)| : x ■ ■ ■ x : Fq 



J DcDc* exp <! Y^ ^^'> 

I nt=0 ns,i 



4(nJ,ai(ns) 
{ns,nt) - Ci{ns,nt + 1)] 



Lt-l 



(21) 



nt=Q 
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where Ci{ns, Lt) = —Ci{ns,0). We use aj(ris) and aj(ns) to denote the fermion annihilation 
and creation operators respectively for spin i. The : symbols in the top line of (l2Tl) indicate 
normal ordering. Let us define the transfer matrix operator 



M„,(s) =: exp 



aUns)ai{ns 



We can write the partition function as 

Zoc j DsTr{ML,-^{s) x ■■■ x Mo(s)}exp |i ^ 



(22) 



(23) 



We now remove the auxiliary field. Integrating over s{n) in the path integral we obtain 



23| 



Z(x j DcDc* exp [-S (c, c*)] 



(24) 



S{c,c*) [c*(n)Q(n + 0) - {I - Qh)c*{n)c,{n)\ 

\c*{n)ci{n + Is) + c*{n)ci{n - Is) 



(e-^--l)(l-6/.)^5: 



c^[n)c^[n)c^[n)ci[n) 



(25) 



This is probably the simplest formulation for computing Feynman diagrams and the most 
convenient for setting the lattice interaction coefficient C. As shown in {l5] the procedure 
for setting C involves summing the bubble diagrams shown in FIG. [H locating the pole in 
the scattering amplitude, and comparing with Liischer's formula for energy levels in a finite 
periodic cube [241. l25l|. 



Ana. 



scatt 



'^scatt 



(26) 



where Ci = —2.837297, C2 = 6.375183. Taking the limit Oscatt ^ oo we get the interaction 
coefficient for the unitary limit. 

We can use (12T|) again to derive the corresponding transfer matrix formalism without 
auxiliary fields, 

Z (xTr {Ml,-i X ■ ■ ■ X Mq} , (27) 
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* «=o,i,..; 1 « * 



FIG. 1: Two-particle scattering bubble chain diagrams 



where 



exp 



-Cat 



l) (1 - 6/i)2 Y.ns a\{ns)a\{ns)a^{ns)ai{ns 



(28) 



Since there is no time-dependent auxiliary field, M„j is the same for each time step rit. This 
formulation is useful in few-body systems with no sign problem. In |t23i] it was used to 
calculate the binding energy of the SU{A) Wigner-symmetric triton. 



V. TRANSFER MATRIX PROJECTION METHOD 



We use the transfer matrix projection method introduced in Hj]. In this paper we 
consider the values = 1,3,5,7,9,11,13. For each spin we fill the momentum states 
comprising j^E'o'^'^) in the order shown in Table 1. We observe that for each value of A^, 
l^f''^^) has zero total momentum and zero total spin. For N = 7 we have the special case 
where I^I/q^"^) is also invariant under 5*0(3, Z) rotations. 



N 


additional momenta filled 


1 


(0,0,0) 


3 


(f,o,o>,(-f,o,o> 


5 


(o,f,o>,(o,-f,o> 


7 


(o,o,f>,(o,o,-f> 


9 


/2n 2-K rv\ / 2-k 2tt n.\ 
\ L ' L ' ^/ ' \ L ' L ' 


11 


/2tt 2-k n.\ 1 2-K 2-k cA 
\ L ' L ' ' \ L ' L ' 


13 


/n 2-K 2k\ /r, 2k 2-k\ 
L ^ L / A^' L ' L / 



Table 1. Filling sequence of momentum states 
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Using the auxiliary field transfer matrix defined in fl22l) we construct 

\^{t), s) = M^^_,{s) X ■ ■ ■ X Mo(s) \¥D , (29) 
where t = ritaf. This is the analog of 

= e-^*|^f'^^^) (30) 

defined above in the continuum notation. We compute spatial correlations of the difermion 
pair operator ip^^fis) = a|(ns)aj^(ns), 



Jds {^ih),s\ 




J 


(dS (^(ti),s|^(t2), 


s)exp{-i^.[s(n)]^} 



G"r(n„ti,t2) = J . (31) 



The inner products in the numerator and denominator are to be defined shortly. We use 
the superscript 'bare' since matrix elements of the composite operators ip'^ and ^/^^^ diverge 
in the continuum limit. We take care of this renormalization by rescaling the correlation 
function, 

G^2{n,MM) = ^G'Jr(n„ti,t2), (32) 

where 

r= hm Gjr(0,ti,t2). (33) 

M„^(s) consists entirely of single-body operators interacting with the background auxil- 
iary field and has no direct interactions between particles. This may not be obvious from 
the complicated form for M„j(s) in (l22l) . To see this more clearly we pretend for the moment 
that the up-spin particles and down-spin particles are all distinguishable. We label 
the newly distinguishable particles with the label X = 1, 2, ■ ■ ■ , 2N — 1, 2N. This error in 
quantum statistics has no effect on the final answer if the initial and final state wavefunc- 
tions are completely antisymmetric in the up-spin and down-spin variables. Since we have 
exactly one particle of each type X, the normal-ordered operator M„j(s) can be factorized 
as a product of terms of the form 

M„,(.) = nM„^(^), (34) 

X 
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where 

+ ^ 5Z [(^xi^s)ax{n, + Is) + a'xi^s)ax{ns - Is) ■ (35) 



If the particle stays at the same spatial lattice site from time step rit to rit + 1, then the 
corresponding matrix element of M^(s) is 



(l-6/i). (36) 



If the particle hops to a neighboring lattice site from time step rit to rif + 1 then the corre- 
sponding matrix element of M^(s) is h. All other elements of M^(s) are zero. 

We can therefore compute the full 2A^-body matrix element as the square of the deter- 
minant of the single-particle matrix elements, 

{^{h), s \^{h), s) ^ Ml,^,{s) X ■ ■ ■ X Mo(s) 

= [detM(s,ti+t2)]% (37) 

Mij{sM + 12) = {pf\Ml^^{s) X ■ ■ ■ X M^{s) \pf) , (38) 

where i,j go from 1 to and ti + ^2 = -^t^t- The states \pf) are the single-particle 
momentum states for the up spins (or down spins) comprising our Slater determinant initial 
and final state |\I'o'"^). The square of the determinant arises from the fact that we have the 
same transfer matrix elements and the same momentum states \pf) for both up and down 
spins. Since the square of the determinant is nonnegative, there is no sign problem. 
We sample configurations according to the weight 



k n 



s{n)f + 2 In [|det M(s, h + ^2) |] } • (39) 



261]. This involves computing 



The updating procedure is done using hybrid Monte Carlo 
molecular dynamics trajectories for 

^(^'P) = ^E(p(^))' + ^(^)' (40) 

n 

where p{n) is the conjugate momentum for s{n) and 

V{s) = i ^ {s{n))' - 2 In [|det M{s, h + ^2) |] • (41) 



2 

n 
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More details of the updating procedure are given in ll|]. For each configuration we compute 
the observable 





^2t(n,)^2(0)|^(t2),s) 




^(t2),s) 



This can be written in terms of the matrix M{s, ti + as we now show. 
We start with 

= «^^| Ml,_^{s) X ■ ■ ■ X M^,^{s)4j'\ns)^\0)Mr,,^_,{s) x ■ ■ ■ x Mo(s) l^^^-^^) (43) 

where ti +t2 = Lfat and ^2 = ^t^at. For each spin the '(/'^"'"(ns)^/'^(0) operator replaces one 
matrix element from M(s, ti + ^2), call it the entry in the fcth row and /th column, with the 
new matrix element 

gki{n,,s,hM) = {Pk\ Ml_,{s) X ■■■ X M„^^(s)a5,(n,)ax(0)M„^^_i(s) x ■■■ x M^{s) \pf). 

(44) 

Instead of the full N x N matrix determinant detM(s, ti + ^2)? in this case we get 
Okii'^s, s,ti,t2) times the entry in the fcth row and /th column of the cofactor matrix of 
M{s, ti + t2). This cofactor matrix element is (—1)^+' times the determinant of the remain- 
ing (A^ — 1) X {N — 1) matrix with the kth row and /th column of M{s, ti + ^2) deleted, and 
it can be calculated as 

det M(s, h + t2) X [M~\s, ti + t2)] , (45) 

where is the matrix inverse of M. Summing over k and / and squaring the result for 
the two spins, we get 

{^ih),s\^'\n,)^'i0)\^it2),s) 

= ldetM{s,h+t2)xJ29kiins,s,h,t2)[M-\s,h+t2)]A • (46) 

I k,l ) 

Therefore our observable is 

0{nsMM.s) = \^gki{ns,sMM)[M-\sM + t2)\^}^ . (47) 

By measuring the ensemble average of 0{ns,ti,t2, s) we get an unbiased estimate of 
Gjr(ns,ti,t2). 
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VI. LATTICE PARAMETERS AND ERROR ESTIMATES 



Our choice for the physical values of the fermion mass and lattice spacings are irrelevant 
to the universal physics of the unitary limit. Nevertheless we must assign values to these 
parameters, and the values we choose are motived by the dilute neutron system. We use 
a fermion mass of 939 MeV and lattice spacings a = (50 MeV)~^, at = (24 MeV)~^. For 
these parameters we find in the unitary limit, Cphys = —1-203 x 10""^ MeV~^. 

For each simulation we compute roughly 2 x 10^ hybrid Monte Carlo trajectories, split 
across four processors running completely independent trajectories. Averages and errors are 
computed by comparing the results of each processor. We use double precision arithmetic 
to compute det M{s, ti + ^2) and 0{ns, ti,t2,s). All systematic errors produced by double 
precision roundoff error and exceptional configurations are monitored in the following way. 
We introduce a small positive parameter e and reject any hybrid Monte Carlo trajectories 
which generate a configuration with 

|detM(s,ti + t2)| < e"^ n \Mn{s,ti + t2)\. (48) 

1=1,.. .,N 

Wc then take the limit e 0"^ to determine if poorly-conditioned matrices make any de- 
tectable contribution to our observables. We consider values for e as small as 10~^. If as 
we take e — > 0+ any systematic error can be detected above the stochastic error level, then 
we throw out the measurement and do not include it in the final results. The error bars 
we present are therefore estimates of the total error for each lattice system. There are no 
additional errors other than the lattice spacing dependence. 

In the unitary limit our Euclidean variables ti and t2 can be replaced by the dimen- 
sionless combinations and More convenient though is to use the dimensionless 

combinations Epti and Ept2, where Ep is the Fermi energy 

At unitarity we can reach the continuum limit at fixed particle number by increasing L, 
the length of the periodic cube in lattice units. This may seem an unusual way to take 
the continuum limit. It works only in scale invariant theories such as the unitary limit or 
noninteracting fermions where the only physical scale is the interparticle spacing. Since 
the number of particles is fixed, the spacing between particles as measured in lattice units 
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increases as we increase L. Therefore L — oo corresponds with the continuum hmit. By 
comparing results for different L we obtain an error estimate for the extrapolation to the 
continuum limit. 

VII. NUMERICAL CROSSCHECKS 

To test the lattice codes, we have run simulations for several systems where the final 
answer can be calculated accurately by alternative means. For the first test we consider the 
noninteracting fermion system for the 7, 7 system with total momentum P = and total 
spin S = 0. We take L = 4 and set Lt large enough to extract the limit 

lim GlTin,,h,h) (50) 

for fis = {rix, 0, 0) with rix = 0,1,2. We perform the numerical check using the path integral 
action S{c,c*) in (!25|) . For temperature T = 0.022Ep and chemical potential /i = 0.97Ef 
we apply the free field Feynman rules for the ip"^ correlation function. For this chemical 
potential at such low temperatures we should see 7 up spins and 7 down spins in the ground 
state with P = and 5 = 0. A comparison of the simulation results and the free grand 
canonical calculations is shown Table 2. We see that the free fermion results agree to 
six-digit precision. 





Simulation results 


Free grand canonical 





1.19629 X 10-2 


1.19629 X 10-2 


1 


6.10352 X 10^3 


6.10352 X 10-3 


2 


2.19727 X 10-3 


2.19727 X 10-3 



Table 2. Simulation results and free grand canonical calculations for 7, 7 

Next we turn on a weak attractive coupling Cphys = —1.25 x 10"^ MeV-^ for the same 7, 7 
system with P = and S = 0. This coupling corresponds with a weak-coupling low-density 
expansion parameter kpagcatt = —0.043. Using T = 0.022Ep and fi = 0.97Ep we use the 
path integral action S{c,c*) to compute the free fermion result as well as the O (kpascatt) 
correction to the tp"^ correlation function. The O {kpascatt) term requires summing the 
two-particle bubble chain shown in FIG. [U A comparison of the simulation results and 
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perturbative calculations is shown in Table 3. The O [kpa'^^att) correction should be of size 
roughly 10~^ for = 0, 5 x 10~^ for = 1, and 2 x 10~^ for = 2. We see that the 
simulation results and perturbative calculations agree within deviations roughly matching 
the size estimates for the O (^I'dgcatt) term. 



rix 


Simulation results 


Perturbative grand canonical 





1.3715(1) X 10-2 


1.3775 X 10-2 


1 


6.999(2) X 10-3 


7.0295 X 10-3 


2 


2.536(1) X 10-3 


2.5476 X 10-3 



Table 3. Simulation results and perturbative grand canonical calculations for 7, 7 

To test that the code is running properly at the unitary limit, we perform simulations for 
the 1, 1 system at P = and S" = at unitarity. We use L = 4 and compute G^^'^^fis, ^1,^2)- 
We take = {n^, 0, 0) for = 0, 1, 2 and various numbers of time steps rit-^ = rit^. Since 
there are only two particles with zero total momentum there should be no dependence on 
the spatial coordinate n^. For the numerical check we use M„j, the transfer matrix without 
auxiliary fields defined in (1251) . to compute the exact two-body transfer matrix in the rest 
frame. The results for the lattice simulation and the exact two-body calculation are shown 
in Table 4. We see that the simulation results agree with the exact results up to errors the 
size of the estimated stochastic noise. 





nt, = = 2 




rtt^ =nt2 = Q 





6.946(9) X 10"^ 


1.156(4) X 10-3 


1.564(4) X 10-3 


1 


6.88(3) X 10"^ 


1.19(4) X 10-3 


1.58(6) X 10-3 


2 


6.92(1) X 10"^ 


1.147(6) X 10-3 


1.53(3) X 10-3 


Exact 


6.948 X 10"^ 


1.153 X 10-3 


1.573 X 10-3 



Table 4: Simulation results and exact two-body calculation for 1, 1 at unitarity 
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FIG. 2: One insertion of the operator a]{r)al{r)a'^{r)ai{r) 



VIII. RESULTS FOR THE RENORMALIZATION CONSTANT T 

Throughout we consider systems with total momentum P = and total spin S = 0. The 
renormalization constant F for the N, N system is given by 

F = lim Gjr(0,ti,t2) = (^ol a|(0)a{(0)at(0)ai(0) |^o) (51) 

where |\E'o) is the normalized ground state. F gives the probability that a given lattice 
site has both an up-spin and down-spin particle. As shown in FIG. |2l one insertion of the 
operator a|(r)a| (7^)0-1- (r)aj(r*) produces two extra fermion bubbles. In the continuum limit 
each of these bubbles are proportional to the inverse lattice spacing plus momentum- 
dependent terms which are finite as a — > 0. Therefore in the continuum limit Fphys is 
proportional to plus subleading terms proportional to a~^. 

The coefficient of a^^ in Fphys contains some useful information about physics near the 
unitarity point. We note that 

F = (xj/ol a\{0)a\{0)a^{0)a^{0) l^j/o) ^ (*o| M^, , (52) 

where M„j is the transfer matrix without auxiliary fields defined in (l28l) . and 

^,_(e-^-0(i^, (53) 

at 

We have dropped terms which are 0{at) and vanish as the temporal lattice spacing goes to 
zero. If Eq is the ground state energy then 

(^o|M„J*o) = e"^°"*, (54) 
r ~ I 'Le-^oat _ 
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Let us parameterize the ground state energy per particle near the unitary hmit as an ex- 

Eq 



pansion in kp^a, ^ 



scatt 5 



F "scatt 



N + N 5 2m 
The renormahzation condition relating C and agcatt gives [l 

d m d mC^ d 



(56) 



da. 



scatt 



47r dC'-i 



47r dC> 



Using 



1/3 



and combining (l55l) . (156|1 . (157|) we find 

5 m (ii?o 



^1 



5m2C"2L4 



(57) 



(58) 



(59) 



127r(67r2)^/='Ar4/3 

Since C'pj^y^ is proportional to the lattice spacing a, we deduce that the leading divergence 
of Tphys is proportional to a^^ as predicted before. This anomalous dimensional scaling can 
also be seen from the dependence of F rather than the naive scaling expected for 
an operator which is the square of a local density. 



We determine V by fitting G^'^^(0,t/2,t/2) at large Ept to the asymptotic form 



G'lf%6,t/2,t/2) -be 



-Tj-Ept 



(60) 



The value of T is then used to determine ^i. We show the results for = 3, 5, 7, 9, 11, 13 
and L = 4, 5, 6 in Table 5. We have extrapolated linearly in as L — > oo to remove 
the subleading dependence in G^^^{0,t/2,t/2). If there are no significant changes for 
> 13 we estimate that in the limit oo, = 1.0(1). 



L 


3,3 


5,5 


7,7 


9,9 


11, 11 


13,13 


4 


0.696(2) 


0.647(2) 


0.597(2) 


0.595(2) 






5 


0.77(1) 


0.719(5) 


0.662(3) 


0.661(2) 


0.652(2) 


0.639(2) 


6 


0.84(3) 


0.785(4) 


0.69(1) 


0.711(4) 


0.71(1) 


0.70(1) 


oo 


1.08(4) 


1.05(4) 


0.91(4) 


0.93(4) 


1.00(5) 


1.01(5) 



Table 5: Results for 
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10 12 



FIG. 3: A comparison of (0, i/2, t/2) and the asymptotic fit 1 - ^e"''--^^* for the 9,9 system 
with L = 4, 5, 6. 



In FIG. Owe show a comparison of the renormalized correlation function 

1 



^^2(0,^/2,^2) = ^Gba-(0,t/2,t/2) 



and the asymptotic fit 



1 _ _p-ri-EFt 
F 



(61) 



(62) 



for the 9, 9 system and L = 4, 5, 6. In the unitary limit we expect agreement for different 
values of L when plotted as functions of Ept, and this appear to be the case. 

The energy r] ■ Ep characterizing the exponential decay of the transient signal in 
(0, t/2, t/2) is similar to the energy 5 ■ Ep measured in the function C,N,N(t) described in 



ll| . This might perhaps be a threshold for a certain type of excitation or a maximum in 



the overlap of |^'o''^) with the spectral density projection operator. The answer is not clear. 
In the next section we find another excitation with the same quantum numbers but much 
lower energy. This suggests that our fit with only one exponential time constant may not 
be a reliable method to determine the higher excitation energy, and a multistate analysis 
should be used. The interpretation of the apparent energy scale i] ■ Ep will require further 
study. 
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IX. RESULTS FOR THE t-DEPENDENT PROFILE OF G^2 

We study the dependence of G^2 on the spatial separation rig and Euchdean projection 
time t. In all cases we consider systems with total momentum P = and total spin 5 = 0. 
In FIG. m we show the renormalized correlation function G^2[ns,t/2,t/2) for the 13,13 
system at unitarity with L = 6, Ug = {nx,0,0), and Ept between and 6.22. In order to 
make the periodicity of the lattice visually clear we show two full lattice lengths. 

There are several features in this plot which indicate some interesting physics. First 
of all ip"^ appears to have long-range order. With a total of 26 particles we cannot probe 
distances far beyond kp^. But for large Ept we do find that G^2(ns,t/2, t/2) averaged over 
the entire lattice is greater than 0.5. We recall (7^2 (n^, t/2, t/2) is normalized so that the 
peak value at = is 1 as Ept —>■ oo. The second point of interest is that for Ug the 
transient signal has a very long time constant. This long constant time is not apparent at 
fis = 0. There the transient signal has a shorter time constant, [rj ■ Ep)~^ in the notation 
of the previous section. The slow transient signal corresponds with a much lower energy 
scale and appears to have an approximate cos{27mx/L) — 1 spatial dependence. 

The time constant is easier to see if we plot G^2(n^, t/2, t/2) as a function of Ept. In 
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FIG. 5: The renormalized correlation function (n^, t/2, t/2) versus Ept for the 7, 7 system with 
L = 4, 5, 6 and = (n^:, 0, 0). For visual clarity datapoints with errorbars exceeding 0.1 are not 
shown. 

FIG. [5] we show the renormalized correlation function G^2[ns,t/2,t/2) versus Ept for the 
7, 7 system with L = 4, 5, 6 and = {rix, 0, 0). We can see quite clearly the fast exponential 
tail of the transient signal for Ug = and the slow exponential tail for Ug ^ 0. In the unitary 
limit, G^iins^tjl^tjl) as a function of Ept for different L should agree for the same ratio 
fis/L. The data in FIG. [5] shows quite clearly this unitary limit scaling. 

In FIG. [6] the renormalized correlation function G^2[ns,t/2,t/2) is shown for the 5,5 
system with L = 4, 5, 6 and = (n2;,0,0). We again see the fast exponential tail of the 
transient signal at = and the slow exponential tail for 7^ 0. There is good agreement 
for different values of L with the same fis/L. 

In FIG. [7] G^2(ns, t/2, t/2) is shown for the same 5,5 system with L = 4,5,6, but this 
time along the 2;-axis, Ug = {0,0,nz). In this case the long time constant previously 
seen for 7^ is no longer visible. The distinction between the x- and z-axes can be 



<0,0,0>, L = 4 I — I — I 
<L/4,0,0>, <3L/4,0,0>, L = 4 x-j - 
<2L/4,0,0> ,L = 4 : - « 

<0,0,0>, L = 5 : □ : 

<L/5,0,0>, <4L/5,0,0>, Z. = 5 o - 
<2L/5,0,0> ,<3L/5,0,0>, 1 = 5 : ^ : " 
<0,0,0>, L = 6 '- -■ - I 
<L/6,0,0>, <5L/6,0,0>, 1 = 6 ; • 
<2L/6,0,0> ,<4L/6,0,0>, L = 6 : 

<3L/6,0,0>, L = 6 ^-T-J 




J I I L 
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FIG. 6: The renormalized correlation function G^2 (n^, t/2, t/2) versus Ept for the 5, 5 system with 
L = 4, 5,6 and = {nx,0,0). For visual clarity datapoints with errorbars exceeding 0.1 are not 
shown. 

explained by the non-5'0(3, Z) invariant momentum filling for |^q*^^). As shown in Table 
1, 1^0^^) contains the momentum states (^,0,0) and (—^,0,0) but not the momentum 
states (O, 0, ^) and (O, 0, —^)- This produces an overlap with some low-energy excitation 
with an approximate cos(27™^/L) — 1 profile in G^2{ns, t/2, t/2) but a much weaker overlap 
if the same excitation is aligned along the 2;-axis. 

X. RESULTS FOR THE LOWEST EXCITATION ENERGY 

We now measure the energy of the lowest excitation. We first construct combinations of 
(fis, ti , ^2) which maximize the signal to noise ratio. For N — 3 the low-energy excitation 
does not couple strongly to G^2(ns, ti, ^2) for along the y- and z-axes. Therefore we use 



<0,0,0>, L = 4 I — I — I 

<zy4,0,0>, <3L/4,0,0>, L = 4 '— x--^ - 
<2L/4,0,0> ,L = 4 

<0,0,0>, L = 5 □ i 

<L/5,0,0>, <4L/5,0,0>, L = 5 ^ -o- -< 

<2zy5,0,0> ,<3L/5,0,0>, L = 5 ; a - 
<0,0,0>, L = 6 - i 

<zy6,0,0>, <5zy6,0,0>, L = 6 ; • 

<2zy6,0,0> ,<4zy6,0,0>, L = 6 ; ^ : . 

<3zy6,0,0>, Z. = 6 T ^ 
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FIG. 7: The renormalized correlation function G^2 (n^, t/2, t/2) versus Ept for the 5, 5 system with 
L = 4, 5, 6 and = (0, 0, n^). For visual clarity datapoints with errorbars exceeding 0.1 are not 
shown. 



only data along the x-axis. We define 



= L' ^G^2 (0, h, t2) - - [6-^2 ((1, 0, 0) , ti, t2) + ((-1, 0, 0) , ti, t2)] 
For N — 5 the coupling is weak only along the 2;-axis and so we define 



^L^{G^2{0,t,,h)-- 



G^2{{1, 0, 0) , h, t2) + G^2{{-1, 0, 0) , ti, ts) 
+G^2{{0, 1, 0) , h, t2) + G'^2((0, -1, 0) , ti, ^2) 



(63) 



(64) 
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FIG. 8: One insertion of the operator ^a|(r)aj(r)V'^ [flt(^)^J.(^] • 



For N > 7 we define 



(7^2(0,^1,^2) 



1 

6 



G^2((l, 0, 0) , ti, ts) + ^^^((-l, 0, 0) , ti, t2) 
+G^2((0, 1, 0) , ti, ts) + ^^^((O, -1, 0) , ti, ts) 
+6-^2 ((0, 0, 1) , ti, t2) + ((0, 0, -1) , ti, 



(65) 



If no additional ultraviolet renormalization is required, then in the continuum limit we 
have 



5lG^2{ti,t2) oc -9^ ^^2 (0,^1,^2), 
5%G^2{tiM) oc - {dl + d^y) G^2(0,ti,t2), 
5%,G^2{t^,t2) ^-{dl + dl + G^2(0,ti,t2) 



(66) 
(67) 
(68) 



To prove that no additional renormalization is needed it suffices to show that one insertion 
of the operator 

ia|(f)a{(f)V2 h(r)ai(f)] (69) 

is finite. Just as we found for the insertion of a|(r)a|(r)a|(r)aj(r), the factor of F~^ takes 
care of the divergences from the two extra fermion bubbles. In the unitary limit the two- 
particle Green's function has the form 

1 



(70) 



m\ mpQ — 



where po is the total energy and p is the total momentum of the two fermions. Consider now 
a diagram such as the one shown in FIG. [81 Let Pq be the energy and p be the momentum 
of the internal loop. We cutoff the momentum at A = 7ra~^ and cutoff the energy at /m. 
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If we now count powers of A we get from dpod^p, A~^ from the two Green's functions, 
A~^ from the three fermion propagators, and A^ from the V^. The final result is A~^ and 
so the integral is finite. In this way one can show that all diagrams involving one insertion 
of 

^a\{f)a\{r)V'[a,{r-)a,{r)] (71) 

are ultraviolet finite. This means that in the unitary limit each of the functions S'^G^2 {ti,t2), 
6lyG^2(ti,t2), 6ly^G^2(ti,t2) should be independent of L when considered as functions of 
Epti and i?i?t2- 

In FIG. in] we plot In [SlG,^2 (ti, ^2)] for the 3, 3 system with L = 4,5, 6. We have produced 
data for ti = ^2 and data for t2 fixed at Ept2 = 1.2. The agreement for different values 
of L provides a consistency check of unitary limit scaling. In the continuum language we 
are calculating the logarithm of — 9^G^2(0, ti, ^2)- Using the asymptotic expansion (fTTj) we 
have 

-dlG^2{0,t„t2) = -a^^o(O) -9^Aoi(0)e-(^-^»)*^ 

- 92Aio(0)e-(^i-^«)*i - a2An(0)e~(^i-^»)(*i+*2) + ■ • ■ . (72) 

The fixed ^2 data as ti — » 00 is useful in extracting Ei — Eq since the time dependence must 
be proportional to e-i^i-^o)ti pj^g constant. We see from the plot that In [5^6*^2 (ti, t2)] 
is nearly a straight line for large ti. This indicates a small asymptotic value at ti = 00, 
and a nearly pure exponential signal in this time window. We have fitted a straight line 
to determine the slope of In [(5^(7^2 (ti, 1(2)] with respect to Epti for the fixed t2 data. We 
do the fit with a common slope for L = 4,5,6 but possibly different intercepts for the three 
values of L. While we are only fitting the fixed ^2 data, we see quite clearly the same slope 
appears in the ti = ^2 data. 

In FIG. [To] we plot In [6lyG^2{ti, ^2)] for the 5, 5 system with L = 4,5, 6. We show data 
for ti = t2 and ^2 fixed at Ept2 = 1.6. We do the same linear fits for the t2 fixed data in 
order to extract Ei — Eq. In FIG. Owe plot In [5ly^G^2(ti,t2)] for the 7, 7 system with 
L = 4,5, 6. In FIG. [12] we plot In [Sly,G^2{ti,t2)] for the 13, 13 system with L = 5,6. 

In all cases we find agreement for different values of L as predicted by unitary limit 
scaling. In all cases we also find agreement between the slope of the fixed t2 data and the 
ti = t2 data when plotted as functions of Ep{ti + 12). This implies that the ^-{^^-^0)^1 g^^fj 
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L = 4, f2 fixed 
L = 5, ?2 fixed 
L = 6, ?2 fixed 
L = 4 fit 
L = 5 fit 
L = 6 fit 




FIG. 9: Plot of In [5lG^2{ti,t2)\ for the 3,3 system with L = 4,5,6. We show data for ti = t2 
and t2 fixed at EFt2 = 1.2. 

^-(Ei-Eo)t2 f^QYms are small so that 



-V'G^2{0,ti,t2) ^ -V%o(0) - V%i(0)e 



-iEi-Eo)iti+t2) 



(73) 



with replaced by for the 3, 3 system and d"^ + dy for the 5, 5 system. 

For the moment let us assume that the excitation at energy Ei can be described as 
two unknown constituents moving in opposite directions with momentum the minimum 
nonzero momentum possible on the lattice. This interpretation of the excitation is probably 
an oversimplification, but it serves as a reasonable starting point to compare with known 
excitations such as pairs of phonons, quasiparticles, or rotons. In FIG. [13] we plot [Ei — 
Eq)/ Ep versus momentum k/kp of the unknown constituents for = 3, 5, 7, 9, 11, 13. For 
comparison we show the linear part of the two-phonon dispersion relation at unitarity. We 
use the result [27 1 

for the speed of sound and use the value ^ = 0.25(3) reported in [111]. We see that the 
minimum in {Ei — Eq)/Ef at k ^ 0.8kp falls well below the linear extrapolation for two 
phonons at A; ~ 0.8kp. This appears to rule out a two phonon interpretation. 

Since O.Skp is close to kp, another reasonable interpretation is that the excitation consists 
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FIG. 10: Plot of In [6'^yG^2{ti,t2)] for the 5,5 system with L = 4,5,6. We show data for ti = t2 
and t2 fixed at Ept2 = 1.6. 



of two fermionic quasiparticles. However here we encounter a similar problem. The energy 
El — Eq is far below current estimates of 2 A, where A is the even-odd staggering in the 
ground state energy. Fixed-node Green's function Monte Carlo simulations get a value 
A = O.STfSVE^ for 12 - 20 particles in a periodic cube 3 and A = 0.50(3)£:i. for 54 - 66 



particles 



291]. Our own lattice simulations for A are in progress, but there is no evidence 



for A being nearly an order of magnitude lower than the fixed-node Monte Carlo results. 
It could be that the quasiparticles are interacting strongly, but some mechanism would be 
needed to explain the significant lowering in energy. 

A third possibility we consider is that the excitation is a pair of rotons moving in opposite 
directions. The shape of the dispersion curve in Fig. [13] is similar to the phonon-roton 



spectrum in superfiuid ^He. Landau 
in superfiuid ^He, 

m 



30 



3l| first predicted the existence of a roton minimum 



2/i 



R 



(75) 



and Feynman gave a quantum-mechanical description of a roton as an excitation whose 



wavelength is resonant with the local spatial correlations of identical Bose particles [32|, 



33, 



34l |. We briefly summarize the argument below. 

Let 00 (^1) ■ ■ ■ i^n) be the ground state wavefunction for an interacting system of 
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FIG. 11: Plot of In [Sly^G^2{ti,t2)] for the 7, 7 system with L = 4, 5, 6. We show data for ti = t2 
and t2 fixed at Ept2 = 2.0. 



identical bosons with energy Eg. We construct a trial state (fi, • • • , fV) with momentum 

— * 

k defined as 



N 



V'fe (n, • • • , rjv) = J] ^^''"'^ ^ '^o (ri, • • • , fiv) . 



(76) 



i=l 



The static structure factor S{k) for the ground state can be written in terms of the square 
of the norm of ijj^, 

N 



Since 



we get 



TV X S{k) = j d^^f 

l,m=L 

N 



2m' 



(77) 



(78) 



(79) 
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FIG. 12: Plot of In [6'^y^G^2 {ti,t2)] for the 13, 13 system with L = 5,6. We show data for ti = t2 
and t2 fixed at EFt2 = 3.0. 



By the variational principle, the energy E{k) for the lowest excitation with momentum k 
satisfies the upper bound 

This upper bound applies to phonons at small k as well rotons at larger k. For a relatively 
dense system such as ^He a maximum in S{k) occurs near k ^ 271/(1, where d is the average 
spacing between bosons. Therefore a minimum in the upper bound flHOj) occurs at ~ 27r/(i, 
and this may also be reflected in E{k). On the other hand for very dilute Bose systems local 
spatial correlations are weaker and become significant only at distance scales comparable 



to the scattering 



V(™scatt) 



35|. 



ength. In this case the maximum in S{k) has been shown to occur at 



The simple estimate kji ^ 2TT/d works rather well for superfluid He. The particle density 



of superfluid "^He at T = 1.25 K and P = 1.0 atm has been measured to be [36| 



N -3 

— ^ 2.20 X lO^^A . (81) 



Making a rough estimate for the average spacing 



d-^iy] , (82) 
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FIG. 13: Plot of {E\ — Eq)/Ef versus momentum k/kp of the unknown consituents. For 
comparison we show the hnear part of the two-phonon dispersion relation at unitarity. 



we get the prediction 



97r 1 



(83) 



which is close to the direct experimental measurement of the roton minimum at T = 1.26 K 



and P = 1.00 atm 



37|, 



k'j^^ = 1.902A . 



(84) 



We return now to our fermionic spin- 1/2 system with up spins and down spins at 
unitarity. Let d be the average distance between neighboring pairs of particles. We use 
the same approximation, 

rf-(f)''^' = iV-V3L, (85) 

and find that 

2n finV^^ 



^^"aFT75Z=^yJ A:^-1.61fc^. (86) 
Unfortunately this does not describe the minimum of {Ei — Eq)/Ef at /c ~ O.Sfci?. In fact 
the minimum in [Ei — Eq)/Ef is very close to one-half of kn. This would appear to rule 
out the interpretation as two rotons. 
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It is likely that our systems of 6 — 26 particles have too few particles to get an accurate 
reading for the excitation spectrum. Indeed there are not enough particles to probe the 
small k/kp behavior where the lowest excitations are described by phonons. Perhaps this 
low-energy excitation changes character and goes up in energy as we include more particles 
and long wavelength phonons emerge. This is possible. However we should still answer the 
question of what is happening in the 6 — 26 particle system and how it is able to produce 
an excitation this low in energy. 

XI. TWO-PARTICLE CORRELATIONS 

In order to study the properties of the ground state and the unknown excitation in 
further detail we compute two-particle correlations. For the same A^, system with total 
momentum P = and total spin S* = 0, we define the opposite-spin two-particle correlation 
function, 

p,,{rM,h)=Pu{rMM) = (vl/(tO \m{t2)) ' ^^^^ 

and the same-spin two-particle correlation function, 

(^(ti)|:al(r>,(r>{(0)a,(0):|vl/(t,)) 
p,,{rMM)=Pn{rMM) = (vi/(tO \m{t2)) • ^^^^ 

In the auxiliary field transfer matrix formalism we compute these correlation functions using 



and 



P[i{ns,ti,t2) 

I Ds (vl>(ti), s l^ih), s) exp {-I [s{n)f} 

(90) 

Rather than calculating matrix elements of the two-particle operators in (1891) and (!90ll 
directly, we define 

Mie^,ei, 5i) =: exp e^a\{ns)a^{ns) + e|aj(nja^(n5) + (5^a}(0)a|(0) : . (91) 

30 



J Ds (^(ti 


),s\ a\{ns)a^ 


[ns)a\{0)a^{6)\^{h),s)exp{-lY.n 


in)?} 




1 Ds i^ih) 


,s\^ih),s)exp{-lJ:.[sin)f} 


(89) 


J Ds (vl/(ti), 


s\ : a|(ns)aj^ 


[ns)a\{0)a^{0) : \^{t2), s) exp {-'^^2^ 


s{n)f} 




/ Ds {^(t,) 


,s\^{t2),s)exp{-lZnHn)f} 





We extract the operators needed for the two-particle correlations by taking derivatives, 

a{(n,)at(n,)aj(0)ai(0) = lim M{e^,e^,SO, (92) 

: al{ns)a^{ns)al{0)a^{0) := lim M{e^,e^,S^). (93) 

This construction is useful because M(e-f , ej^, 5J itself looks like a transfer matrix with only 
single-nucleon operators. 

Let us define the new single-particle matrix elements with one insertion of M^(e|, Cj, 5j^), 

Mij{s,ti,t2,e^,ei,Si) 

= (pf I m4+„^^_,(s) X . • . X M^^(s)M^(eT, e^, 5^)M^^^_,(s) x • • • x M^{s) \pf) , (94) 
We let 

/ Ds exp {-i J2n [«(^)]^} det M(s, ti,t2, e^, e^, (^J 

p{ti, t2, e^, e^, 5i) = J . (95) 

/ Ds exp{-l^.[s{n)f}detM{s,ti + t2) 

Then 

d d 

Pl^(n„ti,t2) = lim — — p(ti,t2,et,e^,5^). (96) 
Pu{ns:ti,t2) ^ lim ——p{ti,t2,e^,€i,di), (97) 

XII. RESULTS FOR AND 

We compute p^{ns,t/2,t/2) and pii{ns,t/2,t/2) for = 5 and L = 4 using the same 
initial state |^o^'') used previously to calculate G^^^{ns, t/2, t/2). We recall that our initial 
state for = 5 corresponds with single particle momenta (0,0,0), (^,0,0), (— ^,0,0), 
(O, X'O)' (O' "T"'*^)- lattice calculation we have performed several consis- 

tency checks for p^i{ns-i t/2, t/2) and Puiusi t/2, t/2). Summing p^i{ns-i t/2, t/2) over and 
multiplying by Li^ counts the number of ways to select two particles with opposite spins, 

L^Y^PuinsMM - N\ (98) 

The same procedure for p^i{ns, t/2, t/2) counts the number of ways to select two down-spin 
particles without replacement, 

^'E'"u(^.>V2,t/2) = N{N-1). (99) 
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One last consistency check uses the fact that 



a|(0)at(0)a|(0)ai(0) = a|(0)a|(0)at(0)ai(0). 



(100) 



From this we deduce that 



Pn(0,t/2,t/2) = G^r(0,t/2,t/2). 



(101) 



The resuhs of the three consistency checks are shown in Table 6. We see that all ratios are 
in agreement with the required value of 1. 



Ept 


^j:Pnins,t/2,t/2) 


lv(fe)EPu(^s,t/2,t/2) 


Pn(o.*/2,t/2) 





1.00 


1.00 


1.00 


1.23 


1.00(1) 


1.00(1) 


1.00(2) 


2.46 


1.00(1) 


1.00(1) 


1.01(1) 


3.70 


1.00(1) 


1.00(1) 


1.01(1) 


4.93 


0.99(1) 


0.99(1) 


0.99(2) 


6.16 


1.01(1) 


1.01(1) 


1.02(1) 


7.39 


0.98(2) 


0.98(2) 


0.98(3) 



Table 6: Consistency checks for p^^ and p^. Each entry should equal 1. 

We show results for Gj,^ "(n„ t/2, t/2), p||(n„ t/2, t/2), and pu(ra„ t/2, t/2) with r?, in 
the xy-plane in FIG. O The same results for in the xz-plane are shown FIG. [T51 In our 
contour plots the maximum brightness for G^^'^ corresponds with 0.05, while the maximum 
brightness for p^j and p^ corresponds with 0.01. For = 5 the ground state of the free 
Fermi system is degenerate. We have chosen the initial state j^I'o'^'^) with all momenta 
orthogonal to the 2;-direction. This explains the lack of cubic 5*0(3, Z) invariance in the 
x2-plane data, especially for smaller values of Ept. At the largest value Ept = 7.39 the 
xz-plane data shows signs of cubic rotational invariance and is beginning to resemble the 
xy-plane data at Ept = 7.39. 

We see that already for Ept > 1.23 a sharp maximum in appears at Ug = 0. This 
indicates significant spatial correlation between up-spin and down-spin particles due to the 
strong zero-range attractive force. Another indication of this spatial correlation is that for 
larger values of Ept and Hs away from the origin, the shape and strength of p-^ and pn 
are nearly identical. Another visible pattern is that regions of maximum strength for G^^^ 
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FIG. 14: G'^T{ns,t/2,t/2), p||(n^, t/2, t/2), ^^(n,, i/2, t/2) with n, in the xy-plane for iV = 5 



and L = 4. 
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FIG. 15: G'^T{ns,t/2,t/2), p||(n^, t/2, t/2), pii{ns,t/2,t/2) with in the xz-plane for iV = 5 
and L = 4. 



correlate with regions of minimum strength for pjj. Conversely the minima for G^^"^ are 
maxima in p^. This is most likely the direct result of Pauli blocking. 

The crosslike pattern that emerges in G^^^ and pn is the same 27r/L cosine wave excita- 
tion we have been investigating. We note that the same pattern can be reproduced rather 
simply with two up-spin and two down-spin particles. Let l^^^) be the normalized four 
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particle state, 




{-Px)aUpx) + a|(Px)a|(-Pz) a|(0)a|(0) |0) , 



(102) 



Px = 



{ 




(103) 



Similarly let |\l'z) be 



1 



at 



{-Pz)aUp-) + aUp,)dU-p,) a{(0)a{(0) |0) 



(104) 



L<^V2 



Pz = {0,0,— 




(105) 



We define t-independent expressions G^^^{ns, '^x,z) and pj,j,(ris, ^x,^) for the states l^^x) and 



The first column in FIG. [TBI shows G^'^'^ and pjj for I^E'x.). The second column shows G^^'^ 
and pjj^ averaged for and I^I/^). The third column shows the same data as the second 
column but cropped to a region of size d x d, where d ~ N~^^^L is the average spacing 
between particles. For these contour plots the maximum brightness for G^^^ corresponds 
with 0.00125, while the maximum brightness for pn corresponds with 0.001. We note that 
the intensity scale for is a factor of ten lower than that in FIGS. [H] and [T3 This is 
consistent with the sum rule in 0991) since we have reduced A^(A^ — 1) from 20 to 2. The 
ratio of intensities between G^I^ and pn in FIGS. [H] and [15] is a factor of four larger than 
that in FIG. [161 This implies that the actual excitation has stronger spatial correlations 
between up and down spins than our simple plane-wave states {"^x) and 

We find that the d x d region in the third column of f[T6|) describes the transient signal 
in G^^*^ and p^ for large Ept. In fact this holds true for all systems we have checked, 

= 3,5,7,9,11,13. This would suggest that the excitation is some type of quasi-lD 
system of four particles wrapped around the periodic lattice. The d x d window would 
suggest a transverse width of approximately d for the quasi- ID system. Outside the d x d 
region it is more difficult to discern a universal pattern in the data for all A^. 



Gj.r (n„ ^^,^) = (^^,,|al(ns)at(n,)a{(0)a|(0) |^^,^) 
Pll{ns,^x,z) = {^x,z\ ■■ a{(n,)a|(n,)a{(0)a|(0) : \^x,z) ■ 



(106) 



(107) 
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X only X and z x and z 




X 



FIG. 16: The first column shows and pn for \^x)- The second column shows and pn 

averaged for l^'^^) and |^'^). The third column shows the same data as the second column but 
cropped to a region of size d x d. 



XIII. SUMMARY AND DISCUSSION 



We have presented lattice results for spin-1/2 fermions at unitarity, where the effective 
range of the interaction is zero and the scattering length is infinite. We measured the 
spatial coherence of difermion pairs for a system of particles with equal numbers of up and 
down spins in a periodic cube. Using a transfer matrix projection method with auxiliary 
fields, we analyzed both ground state properties and transient behavior due to low-energy 
excitations. We first measured F, the probability that a given lattice site has both an 
up-spin and down-spin particle. From this we were able to deduce that i^i = 1.0(1), where 
,^1 is defined by 

nTJ^ = - ^i^-'-scLt + 0(A:-a- J] . (108) 

We then measured the pair correlation function at nonzero spatial separation and found clear 
evidence of off-diagonal long-range order. This suggests that the ground state is superfiuid 
with s-wave pairing. 

We also found evidence for a low-energy excitation with energy Ei which ranges be- 
tween 0.05Ef and 0.2Ef for = 3,5,7,9,11,13. We note that recent radio frequency 
measurements of the excitation spectra in ^Li at unitarity has found an energy gap in the 
unpolarized system at roughly 0.2Ef [38]. The low value for the excitation energy would 
tend to drive down the superfiuid critical temperature and create a pseudogap phase where 
single fermionic quasiparticles are gapped at A while the new excitation is gapped at a lower 
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FIG. 17: Starting with the ground state, we insert one extra up-spin and one extra down-spin. 



energy scale. An estimate of the effect on the critical temperature requires a better under- 
standing of the energy gap in large systems, the density of states, and possible interactions. 
It is possible though that this excitation could explain the consistently low critical temper- 
ature measured in recent lattice simulations which have specifically looked for off-diagonal 
long-range order 
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4Q|. 



The low-energy excitation corresponding with Ei in systems = 3, 5, 7, 9, 11, 13 appears 
to be inconsistent with an interpretation as a pair of weakly interacting phonons, fermionic 
quasiparticles, or rotons. By examining the pair correlation function as well as two-particle 
density correlations, we find the excitation has the characteristics of a quasi-lD chain con- 
sisting of two up-spin and two down-spin particles aligned along one of the lattice axes. 
We caution that this excitation could be an artifact due to the periodic boundary with a 
relatively small number of particles. Confirmation in lattice systems with more particles 
and different geometries will be needed. The exact mechanism which produces the quasi-lD 
chain is beyond the scope of this study. However we can present here at least one plausible 
explanation. 

Consider the ground state of the unitary limit system with — 1 up spins and A^ — 1 
down spins. Let d be the average separation between particles. We now introduce one 
extra up-spin and one extra down-spin separated by a distance x > d as shown in FIG. [ITl 
Next we evolve the state forward in Euclidean time using the projection operator e~'^*. If 
X ^ ^/t/m ~ d then the projection time is not sufficient to find the true ground state of 
the A^, A^ system. Instead we have two essentially non-interacting quasiparticles, and the 
expectation value of the energy after Euclidean time t will be roughly 2A above the ground 
state energy, where A is the even-odd pairing gap. 
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^ ^c^ 



FIG. 18: The pairs located near a line segment connecting the extra particles favor an orientation 
with the up spins to the left and down spins to the right. 



X = x = L 

I n \ 

V H 

d = m 

FIG. 19: Two up spins and two down spins aligned along one of the lattice axes. 

If however x > \Jt/m ~ rf, then pairs nearby the extra particles may rearrange themselves 
slightly to lower the total energy as shown in FIG. [181 Pairs located near a line segment 
connecting the extra particles favor an orientation with up spins to the left and down spins 
to the right. The tilt angle with respect to the line segment may prefer to alternate from 
one pair to the next. The resulting state could be described as a pair of quasiparticles 
interacting via a loose chain of pairs connecting them. The transverse width of this chain 
would be roughly equal to d. 

Chain configurations without specified endpoints can also be constructed. In FIG. [19] 
we show a four-particle chain extending across the periodic lattice. From the symmetry of 
the four-particle chain configuration we expect the resonance energy to be minimized when 
the average spacing between particles equals L/2. Using a simple estimate for the average 
separation, d ~ L/N^^^, we find that a minimum in energy should occur at ~ 8. This is 
consistent with the A^-dependence of {Ei — Eq)/ Ep shown in FIG. [131 Further studies will 
be needed to see if longer quasi- ID chains exist and if excitations can been found which do 
not wind around the periodic lattice. If longer chains are possible then we might expect 
resonances in a periodic cube for d = L/j for each integer j > 2. This corresponds with 
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^ j^. However the alternation of tilt angle between neighboring pairs will be frustrated 
for odd j 

In our discussion we have tried to put an emphasis on universal observables in the unitary 
limit. These are observables which agree in different physical systems at unitarity, such 
as critical temperature, critical velocity, speed of sound, low-energy excitation energies, etc. 
The reason for the emphasis on universal observables is that different microscopic theories 
can agree on all low-energy physical observables but disagree on wavefunctions and off-shell 
Green's functions. In order to make contact with recent observations of superfluidity with 



trapped ultracold atoms 
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45( 1 ■ it seems important for both theory and 



experiment to identify truly universal properties of the unitary limit. There has been 
some recent work on the non-universality of fermion momentum occupation numbers with 
respect to field redefinitions 
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47( 1 . By analogy similar issues have been raised concerning 



the universality of superfiuid condensate fractions for a low-energy effective theory without 
specification of the fundamental microscopic physics. 
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